#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Wed Aug 18 21:41:40 2021

@author: liqingsimac
"""

'''
###############################################################
###############################################################
import math
import numpy as np
import matplotlib.pyplot as plt
from sympy import *
from sympy.plotting import plot
from sympy.plotting import plot_implicit

points = np.array([[10,20], [20, 75], [45, 15], [70, 53], [90, 30]])
weights = np.array([5, 5.1, 5.2, 4.9, 4.8])

x, y = symbols('x y')

p = None
for i in range(len(points)):
  G = (x-points[i][0])**2 + (y-points[i][1])**2 - weights[i]**2   
  p2 = plot_implicit(G, (x,-50,100), (y,-50,100), show=False, line_color='r')
  if p:
    p.extend(p2)
  else:
    p = p2

p.show()
'''

'''
###############################################################
###############################################################
## Example 2.1.1.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
C=[-1000,-300,-100,-50,-10,-4,-1,-0.2,
   0,0.2,1,4,10,50,100,300,1000]
p = None
for k in range(len(C)):
  G = x**2*y**3-C[k]   
  p2 = syp.plot_implicit(G, (x,-5,5), (y,-5,5), show=False, line_color='b')
  if p:
    p.extend(p2)
  else:
    p = p2
p.show()
p.save('ode-example-2-1-1.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2-1-Green's formula
import numpy as np
import matplotlib.pyplot as plt

fig=plt.figure()
ax=fig.add_subplot(111)
ax.hlines(y=0,xmin=-1.5,xmax=1.5)
ax.vlines(x=0,ymin=-1.5,ymax=1.5)
ax.set_aspect('equal')
t=np.linspace(-np.pi/2,-np.pi*3/2,21)
x=np.cos(t)
y=np.sin(t)
ax.plot(x,y,'-')
ax.text(0.1,-1.1,'A',fontsize=15)
ax.text(0.1,0.9,'B',fontsize=15)
for k in range(21):
    ax.plot([0,x[k]],[0,y[k]],'-o')
ax.text(-0.27,-1.3,'C',fontsize=15)
ax.text(-0.45,-1.25,'D',fontsize=15)
ax.text(0.1,0.1,'O',fontsize=15)
fig.savefig('ode-example-2-1-2.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.2.1.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
C=[0,0.1,1,10,100,-0.1,-1,-10,-100]
p = None
for k in range(len(C)):
  G = y**2-1-C[k]*sy.exp(-x**2)/x/x   
  p2 = syp.plot_implicit(G, (x,-3,3), (y,-3,3), show=False, line_color='b')
  if p:
    p.extend(p2)
  else:
    p = p2
#p.show()
p.save('ode-example-2-2-1.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.2.2.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
C=[-3,-2,-1,0,1,2,3,4]
p = None
for k in range(len(C)):
  G = y**2-(x+C[k])**3   
  p2 = syp.plot_implicit(G, (x,-3,3), (y,-3,3), show=False, line_color='b')
  if p:
    p.extend(p2)
  else:
    p = p2
p.show()
p.save('ode-example-2-2-2.png')
'''

'''
## Example 2.2.3.
import numpy as np
import matplotlib.pyplot as plt

k=1
g=10
m=2
a=np.sqrt(k*g/m)
b=np.sqrt(m*g/k)
v0=[1,2,3,4,5,6,7,8,9]
C=(v0+b)/(v0-b)
t=np.linspace(0,1,21)

fig=plt.figure()
ax=fig.add_subplot(111)

for k in range(len(v0)):
    v=b*(C[k]*np.exp(2*a*t)+1)/(C[k]*np.exp(2*a*t)-1)
    ax.plot(t,v)

ax.hlines(y=0,xmin=0,xmax=1)
ax.vlines(x=0,ymin=0,ymax=10)
ax.set_xlabel('t = time')
ax.set_ylabel('v = velocity')

fig.savefig('ode-example-2-2-3.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.3.1.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
C=[-20,-10,-6,-2,-1,0,1,2,6,10,20]
p = None
for k in range(len(C)):
  G = y-x**4/5-C[k]/x  
  p2 = syp.plot_implicit(G, (x,-3,3), (y,-4,4), show=False, line_color='b')
  if p:
    p.extend(p2)
  else:
    p = p2

p.show()

#p.save('ode-example-2-3-1.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.3.2.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
C=[-100,-10,-1,-0.1,-0.01,-0.001,0,0.001,0.01,0.1,1,10,100]
p = None
for k in range(len(C)):
  G = y + sy.cos(x)/5 - 2*sy.sin(x)/5 - C[k]*sy.exp(-2*x)  
  p2 = syp.plot_implicit(G, (x,-7,7), (y,-3,3), show=False, line_color='b')
  if p:
    p.extend(p2)
  else:
    p = p2
p.show()
p.save('ode-example-2-3-2.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.4.2.
import sympy as sy
#import sympy.plotting as syp
from sympy.solvers import dsolve 
x = sy.symbols('x')
v = sy.symbols('v', cls=sy.Function)
ode1=v(x).diff(x,1)-x*v(x)-sy.sin(x)
sol1=dsolve(ode1,v(x))
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Stewart section 5.5 Example
import numpy as np
import matplotlib.pyplot as plt 
theta=np.linspace(0,2*np.pi,201) 
r1=np.abs(np.cos(5.0*theta) - 1.5*np.sin(3.0*theta)) 
r2=theta/np.pi
r3=2.25*np.ones_like(theta)
fig = plt.figure()
ax = fig.add_subplot(111,projection='polar')
ax.plot(theta,r1,label='trig')
ax.plot(5*theta,r2,label='spiral')
ax.plot(theta, r3,label='circle')
ax.legend(loc='best')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2.4.3
import numpy as np
import matplotlib.pyplot as plt

theta=np.linspace(-2*np.pi,2*np.pi,1001)
r1=np.exp(theta)/100
r2=2*np.exp(theta)/100
r3=3*np.exp(theta)/100
r4=theta
r5=2*theta
r6=3*theta

fig = plt.figure()
ax = fig.add_subplot(121,projection='polar')
ax.plot(theta,r1)
ax.plot(theta,r2)
ax.plot(theta,r3)
bx = fig.add_subplot(122,projection='polar')
bx.plot(theta,r4)
bx.plot(theta,r5)
bx.plot(theta,r6)
plt.tight_layout(pad=0.4, w_pad=0.5, h_pad=1.0)
plt.savefig('ode-example-2-4-3.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2-5-3, d(tan(theta)) vs. d(theta)
import numpy as np
import matplotlib.pyplot as plt

#u=tan(theta)
#theta=arctan(u)

# =============================================================================
# fig = plt.figure()
# ax1 = fig.add_subplot(111)
# theta=np.linspace(-np.pi/2+0.1, np.pi/2-0.1, 22)
# u=np.tan(theta)
# ax1.plot(theta,u,'-')
# afew=[0,3,6,9,12,15,18,21]
# ax1.plot(theta[afew],np.zeros_like(theta[afew]),'.-')
# ax1.plot(np.zeros_like(u[afew]),u[afew],'.-')
# ax1.plot(theta[afew],u[afew],'.')
# ax1.text(0.6,-1,r'$d\theta$',fontsize=12) # t=theta
# ax1.text(-0.5,1,'du',fontsize=12)
# ax1.text(-1,9,r'$u=\tan(t)$',fontsize=12)
# =============================================================================

fig = plt.figure()
ax2 = fig.add_subplot(111)
u=np.linspace(-10, 10, 51)
theta=np.arctan(u)
ax2.plot(u,theta,'-')
afew=np.arange(0,51,5)
ax2.plot(u[afew],np.zeros_like(u[afew]),'.-')
ax2.plot(np.zeros_like(theta[afew]),theta[afew],'.-')
ax2.plot(u[afew],theta[afew],'.')
ax2.text(2.5,0.1,r'$du$',fontsize=12) # t=theta
ax2.text(-1.5,1.15,r'$d\theta$',fontsize=12)
ax2.text(5,1,r'$\theta=\arctan(u)$',fontsize=12)
ax2.text(4.4,0.7,r'$d\theta=\frac{du}{1+u^2}$',fontsize=12)

ax2.plot(u[[30,35]],np.zeros_like(u[[30,35]]),linewidth=3)
ax2.plot(np.zeros_like(theta[[30,35]]),theta[[30,35]],linewidth=3)

fig.savefig('ode-example-2-5-3.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2-6-1
import numpy as np
import matplotlib.pyplot as plt

K=np.arange(start=-2,stop=2.5,step=0.5)
C=np.arange(start=0,stop=5,step=1)

fig = plt.figure()
ax = fig.add_subplot(111)

x=np.linspace(-5,5,51) 
for m in range(len(K)):
    y=K[m]*x**2
    ax.plot(x,y,'b-')

t=np.linspace(0,2*np.pi,51)
for k in range(len(C)):
    x2=C[k]*np.cos(t)
    y2=C[k]*np.sin(t)/np.sqrt(2)
    ax.plot(x2,y2,'r-')

ax.set_xlim([-5,5])
ax.set_ylim([-4,4])

ax.set_aspect('equal')
plt.savefig('ode-example-2-6-1.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2-6-2
import numpy as np
import matplotlib.pyplot as plt

N0=9.7092
k=0.0145
t0=1979

a=0.029
b=a/20
t=np.linspace(1979,2049,18)
N1=N0*np.exp(k*(t-t0))
N2=a*N0*np.exp(a*(t-t0))/(a-b*N0+b*N0*np.exp(a*(t-t0)))

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax1.plot(t,N1,'bo')
ax2 = fig.add_subplot(122)
ax2.plot(t,N2,'bo')

#ax1.set_xlim([0,10])
#ax1.set_ylim([0,20])

plt.savefig('ode-example-2-6-2.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
##Example 2-6-2-plus
import numpy as np
import matplotlib.pyplot as plt

t=np.linspace(-10,10,101)
N1=np.exp(0.5*t)
N2=6/(1+np.exp(-0.5*t))

fig = plt.figure()
ax1 = fig.add_subplot(121)
ax1.plot(t,N1,'b-')
ax2 = fig.add_subplot(122)
ax2.plot(t,N2,'b-')

#ax1.set_xlim([0,10])
#ax1.set_ylim([0,20])

plt.savefig('ode-example-2-6-2-plus.png')
###############################################################
###############################################################
'''

'''
###############################################################
###############################################################
## Example 2.6.3.
import sympy as sy
import sympy.plotting as syp
x, y = sy.symbols('x y')
mu1=1
delta1=0.1
lambda1=1.5
sigma1=0.075
C=2

G=delta1*x + sigma1*y - mu1*sy.log(x) - lambda1*sy.log(y) - C
p=syp.plot_implicit(G, (x,0,100), (y,0,200), line_color='b')

p.save('ode-example-2-6-3.png')
###############################################################
###############################################################
'''

